import os
import sys
module_path = sys.path.append('/home/dpirvu/.local/lib/python3.8/site-packages/')
module_path = sys.path.append('/home/dpirvu/.local/bin/')
for i in sys.path:
i = i.replace("8", "7")
if i not in sys.path:
sys.path.append(i)
print(i)
import numpy as np
import random
import math
import statistics as stat
import matplotlib
import matplotlib.pyplot as plt
import scipy as scp
import scipy.optimize as sco
import scipy.signal as scs
import scipy.special as ssp
import scipy.integrate as sci
from itertools import zip_longest, groupby, cycle
from scipy.ndimage import gaussian_filter1d
from uncertainties import ufloat
np_load_old = np.load
np.load = lambda *a,**k: np_load_old(*a, allow_pickle=True, **k)
nLat = 1024
nTime = nLat*2
nSimTime = nLat*2
lSim = 0
nSims = 100
phi0 = 1.
m2eff = 1.
lamb = 1.
lenLat = 100.;
phi_initial = 0
nCols = 2
nyq = nLat//2
spec = nyq
alpha = 8
step = 1
dx = lenLat/nLat
dk = 2.*np.pi/lenLat
dt = dx/alpha
dtout = dt*alpha*step
plots_file = '/home/dpirvu/project/thermal_bubbles/plots/'
pickle_file = '/home/dpirvu/project/pickle_location/thermal_bubbles/'
sim_location = lambda nT, nL, tem, phi, lam, m2, gam, sim: '/gpfs/dpirvu/thermal_bubbles/phi4_classical_solution_with_fluctuations_m2{:.4f}'.format(m2)+'_gam{:.4f}'.format(gam)+'_t'+str(nT)+'_x'+str(nL)+'_temp{:.4f}'.format(tem)+'_phi0{:.4f}'.format(phi)+'_lamb{:.4f}'.format(lam)+'_sim'+str(sim)+'_fields.dat'
classical_field = lambda gam, tem: np.asarray(np.reshape(np.genfromtxt('/gpfs/dpirvu/thermal_bubbles/phi4_m21.0000_gam{:.4f}'.format(gam)+'_t1024_x1024_temp{:.4f}'.format(tem)+'_phi01.0000_lamb1.0000_fields.dat'), (nTime, nLat)))
smoothen = lambda slice, sigma: np.fft.ifft(np.fft.fft(slice)*Filter(sigma))
def extract_sim(nT, nL, tem, phi, lam, m2, gam, sim, sigma=0.):
data = np.genfromtxt(sim_location((nSimTime if gam >= 0. else nTime), nL, tem, phi, lam, m2, gam, sim))
simulation = [np.reshape(data[:,col], (nT, nL)) if sigma == 0. else np.asarray([smoothen(slic, sigma) for slic in simulation[col]]) for col in range(nCols)]
return np.asarray(simulation)
# field spectrum
# prediction from lattice renormalization
norm = lambda phi: 1./ phi / np.sqrt(2.*lenLat)
freqs = np.fft.fftshift(np.fft.fftfreq(nLat, d=dx))*2*np.pi
omega = np.sqrt(freqs**2. + m2eff)
Filter = lambda sigma: np.exp(-0.5* (freqs*dx*sigma)**2.)
spectral_field = lambda phi, tem: np.asarray([norm(phi) * np.sqrt(2./(np.exp(omega[k]/tem) - 1.)) / np.sqrt(omega[k]) if k != nLat//2 else 0 for k in range(nLat)])
pspec = lambda phi, tem, sigma: np.abs(spectral_field(phi, tem) * Filter(sigma))**2.
sigma0sq_th = lambda phi, tem, sigma: np.sum(pspec(phi, tem, sigma))
V = lambda x, lam, gam: (- 0.5 * m2eff * x**2. + lam * x**4. / 24. - gam * x) if gam == 0.25 else (- 0.5 * m2eff * x**2. + lam * x**4. / 24. - gam * x**3. / 6. if x < 0. else - 0.5 * m2eff * x**2. + lam * x**4. / 24.) if gam == 0.5 else (- 0.5 * m2eff * x**2. + lam * x**4. / 24.) if gam == 0. else (0.5 * m2eff * x**2. + lam * x**4. / 24.)
VV = lambda x: - m2eff + lamb * x**2. / 2.
VVVV = lambda lam: lam
th_masses = lambda phi, lam, gam, tem, sigma: (m2eff if gam == -1. else -m2eff) + 0.25 * VVVV(lam) * sigma0sq_th(phi, tem, sigma)
masses_from_classical_data_ps = lambda lam, gam, tem: (m2eff if gam == -1. else -m2eff) + 0.25 * VVVV(lam) * np.mean(classical_field(gam, tem)[0]**2.)
egtemps = np.linspace(0, max(freqs), 100)
plt.plot(egtemps, np.abs([th_masses(phi0, lamb, 0., tem, 0.) for tem in egtemps]))
plt.xlabel('Maximum Temperature Range')
plt.ylabel('Expected Field Mass')
plt.show()
#gamma = list_gamma[0]
#for temp in list_temp[::100]:
# print('temp, gamma:', temp, gamma)
# field = classical_field(gamma, temp)
# field = np.asarray([np.roll(slice, nLat//2) for slice in field])
# fftfield = np.abs(np.fft.fftshift(np.fft.fft2(field)))**2.
# fig, ax = plt.subplots(1, 2, figsize = (10, 4))
# im = ax[0].imshow(field, aspect='auto', interpolation='none', origin='lower')
# clb = plt.colorbar(im, ax = ax[0])
# im = ax[1].imshow(fftfield, aspect='auto', interpolation='none', origin='lower')
# clb = plt.colorbar(im, ax = ax[1])
# plt.suptitle('Temperature = '+str(temp))
# plt.show()
#gamma = 0
#for temp in list_temp:
# fig, ax = plt.subplots(1, 1, figsize = (20,4))
# plt.plot(freqs, pspec(phi0, temp, 0.), label='theoretical')
# fld = classical_field(gamma, temp)[0]
# plt.plot(freqs, np.abs(np.fft.fftshift(np.fft.fft(fld)))**2./nLat**2, 'o', ms=1, label='data')
# plt.legend(); plt.show()
xlist = np.linspace(-4, 4, 100)
plt.figure()
[plt.plot(xlist, [V(x, lamb, gam) for x in xlist], label=f'V($\phi; \gamma = $%5.2f)'%(gam)) for gam in [-1, 0]]#, 0.25, 0.5]]
plt.legend(); plt.show()
list_gamma = [0.]
list_temp = np.arange(1, 11, 1)/10.
fine_temp_list = np.linspace(min(list_temp), max(list_temp), 100)
filter_size = 0.
all_data = np.asarray([[[extract_sim(nTime, nLat, temp, phi0, lamb, m2eff, gamma, sim, filter_size) for sim in range(lSim, nSims)] for gamma in list_gamma] for temp in list_temp])
#dataPS = lambda tem, sim, time: np.abs(np.fft.fftshift(np.fft.fft(all_data[tem][0][sim][0][time])))**2./nLat**2.
#for tem in range(len(list_temp)):
# time = 0
# fig, ax = plt.subplots(1, 1, figsize = (20,4))
# plt.plot(freqs, pspec(phi0, list_temp[tem], 0.))
# plt.plot(freqs, np.mean([dataPS(tem, sim, time) for sim in range(nSims-lSim)], axis=0), 'o', linestyle='None', ms=2)
## [plt.plot(freqs, dataPS(tem, sim, time), 'o', ms=2) for sim in range(nSims-lSim)]
# plt.show()
sim = random.randint(0, nSims-lSim-1); print(sim)
fig, ax = plt.subplots(len(list_gamma)+1, len(list_temp), figsize = (4*len(list_temp)+4, 4*len(list_gamma)+4))
for gamma in range(len(list_gamma)):
for temp in range(len(list_temp)):
im0 = ax[gamma][temp].imshow(all_data[temp][gamma][sim][0], aspect='auto', interpolation='none', origin='lower')
clb = plt.colorbar(im0, ax = ax[gamma][temp])
ax[gamma][temp].set_title('T = '+str(list_temp[temp])+f', $\gamma$ = '+str(list_gamma[gamma]))
plt.suptitle('Simulation '+str(sim))
plt.show()
temp = len(list_temp)//2
gamma = 0
div = 10
fig, ax = plt.subplots(div, int(np.ceil((nSims-lSim)/div)), figsize = (4*(nSims-lSim)+4, div*4))
io = -1
for sim in range(nSims-lSim):
if sim % div == 0: io += 1
im0 = ax[sim%div][io].imshow(all_data[temp][gamma][sim][0], aspect='auto', interpolation='none', origin='lower')
clb = plt.colorbar(im0, ax = ax[sim%div][io])
ax[sim%div][io].set_title('sim = '+str(sim))
plt.suptitle('T = '+str(list_temp[temp])+f', $\gamma$ = '+str(list_gamma[gamma]), weight=10, size='xx-large', y=1)
plt.tight_layout(); plt.show()